Identification pipeline

To make the signature identification procedure more clear, we drew a flowchart. We name this pipeline “Best Practice” and use its short name BP (or bp).

TODO: add flowchart

Apply to PCAWG catalog matrix

PBS file bp_pcawg.pbs content:

#PBS -l walltime=500:00:00
#PBS -l nodes=1:ppn=50
#PBS -S /bin/bash
#PBS -l mem=300gb
#PBS -j oe
#PBS -M w_shixiang@163.com
#PBS -q slst_pub

# Please set PBS arguments above
cd /public/home/wangshx/wangshx/PCAWG-TCGA
module load apps/R/3.6.1

# Following are commands/scripts you want to run
Rscript call_pcawg_bp.R

R script call_pcawg_bp.R content:

library(sigminer)

tally_X <- readRDS("pcawg_cn_tally_X.rds")
tally_X_noLOH <- readRDS("pcawg_cn_tally_X_noLOH.rds")

sigs <- bp_extract_signatures(
  tally_X$nmf_matrix,
  range = 2:30,
  n_bootstrap = 10,
  n_nmf_run = 10,
  cores = 50,
  cache_dir = "pcawg_bp_pynmf",
  keep_cache = TRUE,
  cores_solution = 30,
  pynmf = TRUE,
  use_conda = TRUE
)

saveRDS(sigs, file = "pcawg_cn_sigs_CN176_BP.rds")

Signature number determination

The code above call 100 NMF runs (10 bootstrap catalogs and 10 NMF runs for each bootstrap catalog) for each signature number. Only 1000 NMF runs (20 bootstrap catalogs and 50 NMF runs for each bootstrap catalog) was applied for 176-component classification approach, which has the same arguments proposed by Degasperi et al. (2020).

First, load the results.

library(sigminer)
solution1000 <- readRDS("../BP/BP_PCAWG_1000_Extraction_Result.rds")
solution100 <- readRDS("../data/pcawg_cn_sigs_CN176_BP.rds")

Next we compare some important measures for signature number determination. We can see 100 NMF runs and 1000 NMF runs can achieve very similar results.

bp_show_survey(solution1000, add_score = F, fixed_ratio = F)

bp_show_survey(solution100, add_score = F, fixed_ratio = F)

Basically, the 5 measures can be classified into 3 categories:

1. silhouette and similarity: these two scores indicate the result stability. The signature number with value decreases sharply is preferred.
2. distance and error: these two scores indicate the difference between raw matrix and reconstructed matrix for signatures. The signature number with lower value is better.
3. pos cor: this is constructed by Prof. Liu and Shixiang to see the average positive correlation between signatures. A signature number should be considered if it has lower value. In practice, we find this measure has referential value sometimes.

More see ?bp_show_survey in R console for how they are calculated.

In previous study, Degasperi et al. (2020) used silhouette and error for signature number determination, PCAWG Mutational Signatures Working Group et al. (2020) used similarity and distance for signature number determination. For plots above we can see that they very close.

From plots above, we observe that the stability drops sharply after 11, so here 11 signatures is selected for this PCAWG dataset.

Similarity between two results with different parameter setting

We have two results above, so we can compare the signatures from two different parameter setting by calculating their pairwise cosine similarity.

sim <- get_sig_similarity(solution1000$object$K11, solution100$object$K11)
ComplexHeatmap::pheatmap(sim$similarity, cluster_rows = FALSE, cluster_cols = FALSE, display_numbers = TRUE)

We can know that the results contain same signature profile! This practice indicates that 100 NMF runs for each signature number is enough.

The Signature object with 11 signature profiles are stored for further analysis.

# Use root directory as work directory
saveRDS(solution1000$object$K11, file = "data/pcawg_cn_sigs_CN176_signature.rds")

Attribution pipeline

There are 3 core steps in this pipeline:

1. Step1: We run n = 100 optimizations on bootstrapped data and obtain a distribution of activities for each sample.
2. Step2: We consider the median activities as point estimates.
3. step3: In each sample, we reset a signature activity to 0 if its activity distribution below the threshold of 1% of total mutations in the sample with one-side T test p value >0.05.

This pipeline is similar to previously reported Huang, Wojtowicz, and Przytycka (2018). However, there are some different settings:

1. We use 1% of total mutations as a threshold here instead 5% in Degasperi et al. (2020), because we observed some copy number segment contributions are removed from results with 5% but they can be used for better explanation of sample copy number profile reconstruction.
2. We use T test for P value calculation here instead of empirical P value calculation based on proportion in Degasperi et al. (2020) for keeping more positive events.

TODO: add flowchart

Apply to PCAWG catalog matrix and signatures

We apply the attribution pipeline above.

# Use root directory as work directory
library(sigminer)
tally_X <- readRDS("data/pcawg_cn_tally_X.rds")
pcawg_types <- readRDS("data/pcawg_type_info.rds")

sc <- pcawg_types$cancer_type
names(sc) <- pcawg_types$sample

pcawg_expo <- bp_attribute_activity(
  solution1000$object$K11,
  sample_class = sc,
  nmf_matrix = tally_X$nmf_matrix,
  return_class = "data.table"
)

saveRDS(pcawg_expo, file = "data/pcawg_cn_sigs_CN176_activity.rds")

Now we can know how well each sample catalog profile can be constructed from the signature combination.

pcawg_expo <- readRDS("../data/pcawg_cn_sigs_CN176_activity.rds")
summary(pcawg_expo$similarity)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.4446  0.8431  0.9043  0.8950  0.9649  0.9999 

hist(pcawg_expo$similarity, breaks = 100, xlab = "Reconstructed similarity", main = NA)

Most of samples are well constructed. Next we will use a similarity threshold 0.75 to filter out some samples for removing their effects on the following analysis.

Signature similarity

Here we explore the similarity between copy number signatures identified in PCAWG data. High cosine similarity will cause some issues in signature assignment and activity attribution (Maura et al. 2019).

sim <- get_sig_similarity(pcawg_sigs, pcawg_sigs)

pheatmap::pheatmap(sim$similarity, cluster_cols = F, cluster_rows = F, display_numbers = TRUE)

Most of signature-pairs have similarity around 0.1.

Signature activity distribution

Transform data to long format.

df_abs <- get_sig_exposure(pcawg_sigs) %>%
  tidyr::pivot_longer(cols = starts_with("Sig"), names_to = "sig", values_to = "activity") %>%
  dplyr::mutate(
    activity = log10(activity + 1),
    sig = factor(sig, levels = paste0("Sig", 1:11))
  )

df_rel <- get_sig_exposure(pcawg_sigs, type = "relative", rel_threshold = 0) %>%
  tidyr::pivot_longer(cols = starts_with("Sig"), names_to = "sig", values_to = "activity") %>%
  dplyr::mutate(
    sig = factor(sig, levels = paste0("Sig", 1:11))
  )

Plot with function in sigminer.

show_group_distribution(
  df_abs,
  gvar = "sig",
  dvar = "activity",
  order_by_fun = FALSE,
  g_angle = 90,
  ylab = "log10(activity+1)"
)

show_group_distribution(
  df_rel,
  gvar = "sig",
  dvar = "activity",
  order_by_fun = FALSE,
  g_angle = 90,
  ylab = "relative activity"
)

Copy number profile for representative samples

Here we will draw copy number profile of 2 representative samples for each signature.

get_sample_from_sig <- function(dt, sig) {
  res <- head(dt[order(dt[[sig]], decreasing = TRUE)], 2L)
  res
}

Read CopyNumber object for PCAWG data.

pcawg_cn_obj <- readRDS("../data/pcawg_cn_obj.rds")
samp_summary <- pcawg_cn_obj@summary.per.sample
rel_activity <- get_sig_exposure(pcawg_sigs, type = "relative", rel_threshold = 0)
rel_activity <- rel_activity[, lapply(.SD, function(x) {
  if (is.numeric(x)) round(x, 3) else x
})]

dir.create("../output/enrich_samples/", showWarnings = FALSE)
for (i in paste0("Sig", 1:11)) {
  cat(paste0("Most enriched in ", i, "\n"))
  s <- get_sample_from_sig(rel_activity, i)
  print(s)
  plist <- show_cn_profile(pcawg_cn_obj,
    samples = s$sample,
    show_title = TRUE,
    return_plotlist = TRUE
  )
  plist <- purrr::map2(plist, s$sample, function(p, x) {
    s <- samp_summary[sample == x]
    text <- paste0(
      "n_of_seg:", s$n_of_seg, "\n",
      "n_of_amp:", s$n_of_amp, "\n",
      "n_of_del:", s$n_of_del, "\n"
    )
    p <- p + annotate("text",
      x = Inf, y = Inf, hjust = 1, vjust = 1,
      label = text
    )
    p
  })
  p <- cowplot::plot_grid(plotlist = plist, nrow = 2)
  print(p)
  ggplot2::ggsave(file.path("../output/enrich_samples/", paste0(i, ".pdf")),
    plot = p, width = 12, height = 6
  )
}

Most enriched in Sig1
     sample  Sig1  Sig2  Sig3  Sig4 Sig5 Sig6  Sig7  Sig8  Sig9 Sig10 Sig11
1: SP105253 0.717 0.045 0.000 0.001    0 0.03 0.183 0.000 0.003 0.003 0.018
2:  SP98941 0.682 0.200 0.002 0.000    0 0.00 0.002 0.099 0.000 0.009 0.005

Most enriched in Sig2
    sample  Sig1  Sig2  Sig3  Sig4  Sig5  Sig6  Sig7  Sig8  Sig9 Sig10 Sig11
1: SP55235 0.056 0.796 0.000 0.006 0.007 0.031 0.010 0.066 0.000     0 0.027
2: SP56303 0.019 0.779 0.004 0.000 0.000 0.102 0.064 0.026 0.005     0 0.001

Most enriched in Sig3
     sample Sig1 Sig2  Sig3  Sig4  Sig5 Sig6  Sig7 Sig8  Sig9 Sig10 Sig11
1: SP112911    0    0 0.988 0.002 0.006    0 0.000    0 0.004     0     0
2: SP112817    0    0 0.959 0.006 0.013    0 0.005    0 0.016     0     0

Most enriched in Sig4
     sample  Sig1  Sig2  Sig3  Sig4  Sig5 Sig6  Sig7  Sig8  Sig9 Sig10 Sig11
1: SP114903 0.144 0.000 0.005 0.615 0.187    0 0.008 0.003 0.030 0.008 0.000
2: SP114898 0.166 0.001 0.005 0.594 0.194    0 0.006 0.002 0.022 0.006 0.002

Most enriched in Sig5
     sample Sig1 Sig2 Sig3 Sig4 Sig5 Sig6 Sig7 Sig8 Sig9 Sig10 Sig11
1: SP105708    0    0    0    0    1    0    0    0    0     0     0
2: SP107359    0    0    0    0    1    0    0    0    0     0     0

Most enriched in Sig6
     sample  Sig1  Sig2  Sig3  Sig4  Sig5  Sig6  Sig7  Sig8 Sig9 Sig10 Sig11
1: SP116525 0.033 0.021 0.001 0.002 0.000 0.905 0.004 0.028    0 0.003 0.004
2: SP121824 0.000 0.000 0.000 0.013 0.028 0.800 0.001 0.000    0 0.002 0.156

Most enriched in Sig7
     sample  Sig1  Sig2 Sig3 Sig4 Sig5  Sig6  Sig7  Sig8  Sig9 Sig10 Sig11
1: SP116345 0.006 0.349    0    0    0 0.003 0.591 0.051 0.000 0.000     0
2: SP124197 0.069 0.104    0    0    0 0.000 0.591 0.156 0.028 0.053     0

Most enriched in Sig8
    sample  Sig1 Sig2 Sig3 Sig4 Sig5 Sig6 Sig7  Sig8 Sig9 Sig10 Sig11
1:  SP5381 0.000    0    0    0    0    0    0 1.000    0     0     0
2: SP78664 0.014    0    0    0    0    0    0 0.986    0     0     0

Most enriched in Sig9
     sample Sig1 Sig2  Sig3  Sig4  Sig5 Sig6 Sig7 Sig8  Sig9 Sig10 Sig11
1: SP106797    0    0 0.073 0.029 0.136    0    0    0 0.762     0     0
2: SP135354    0    0 0.043 0.000 0.251    0    0    0 0.705     0     0

Most enriched in Sig10
     sample  Sig1 Sig2 Sig3  Sig4  Sig5 Sig6  Sig7  Sig8  Sig9 Sig10 Sig11
1: SP117086 0.000    0    0 0.000 0.038    0 0.065 0.094 0.043 0.760     0
2: SP117223 0.071    0    0 0.007 0.197    0 0.000 0.015 0.022 0.688     0

Most enriched in Sig11
     sample  Sig1  Sig2  Sig3  Sig4  Sig5  Sig6  Sig7  Sig8  Sig9 Sig10 Sig11
1:  SP11235 0.047 0.030 0.051 0.140 0.046 0.002 0.051 0.028 0.056 0.065 0.485
2: SP112194 0.073 0.128 0.003 0.032 0.031 0.054 0.032 0.091 0.012 0.060 0.484

Ideograph for copy number profile reconstruction

TODO: select one/two samples and draw copy number profile -> catalog profile -> signature relative activity. Similar to Nature paper.

Signature activity in each cancer type

Here we draw distribution of a signature across cancer types.

show_group_distribution(
  df_abs,
  gvar = "cancer_type",
  dvar = "Sig1",
  order_by_fun = FALSE,
  g_angle = 90,
  point_size = 0.3
)

We have many signatures here, so we output them to PDF files.

dir.create("../output/cancer-type-dist", showWarnings = F)
signames <- paste0("Sig", 1:11)
for (i in signames) {
  pxx <- show_group_distribution(df_abs,
    gvar = "cancer_type",
    dvar = i, order_by_fun = FALSE,
    ylab = i,
    g_angle = 90, point_size = 0.3
  )
  ggplot2::ggsave(file.path("../output/cancer-type-dist/", paste0("Absolute_activity_", i, ".pdf")),
    plot = pxx, width = 12, height = 6
  )
  pxx <- show_group_distribution(df_rel,
    gvar = "cancer_type",
    dvar = i, order_by_fun = FALSE,
    ylab = i,
    g_angle = 90, point_size = 0.3
  )
  ggplot2::ggsave(file.path("../output/cancer-type-dist/", paste0("Relative_activity_", i, ".pdf")),
    plot = pxx, width = 12, height = 6
  )
}
rm(pxx)

Signature landscape

Define a signature which is detectable if this signature contribute >5% exposures in a sample.

df <- df_rel %>%
  dplyr::mutate_at(dplyr::vars(dplyr::starts_with("Sig")), ~ ifelse(. > 0.5, 1L, 0L)) %>%
  tidyr::pivot_longer(
    cols = dplyr::starts_with("Sig"),
    names_to = "sig", values_to = "detectable"
  )

df2 <- df_rel %>%
  tidyr::pivot_longer(
    cols = dplyr::starts_with("Sig"),
    names_to = "sig", values_to = "expo"
  )
df <- dplyr::left_join(df, df2,
  by = c("sample", "cancer_type", "sig")
)

df_type <- df %>%
  dplyr::group_by(cancer_type, sig) %>%
  dplyr::summarise(
    freq = sum(detectable), # directly use count
    expo = median(expo[detectable == 1]),
    n = n(),
    label = paste0(unique(cancer_type), " (n=", n, ")"),
    .groups = "drop"
  )

mps <- unique(df_type[, c("cancer_type", "label")])
mpss <- mps$label
names(mpss) <- mps$cancer_type

summary(df_type$freq)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.000   0.000   0.000   3.395   1.000 113.000 

Show copy number signature landscape.

library(cowplot)

p <- ggplot(
  df_type,
  aes(x = cancer_type, y = factor(sig, levels = paste0("Sig", 1:11)))
) +
  geom_point(aes(size = freq, color = expo)) +
  theme_cowplot() +
  ggpubr::rotate_x_text(60) +
  scale_x_discrete(breaks = mps$cancer_type, labels = mps$label) +
  scale_size_continuous(
    limits = c(5, 230),
    breaks = c(5, 20, 50, 100, 200)
  ) +
  scale_color_gradient2(low = "white", mid = "purple", high = "black", midpoint = 0.3) +
  labs(
    x = NULL, y = "Copy number signatures",
    color = "Median exposure\ndue to signature",
    size = "Tumors with\nthe signature"
  )
p

Signature activity association

For absolute signature activity.

show_cor(pcawg_expo$abs_activity[, -1])

For relative signature activity.

show_cor(pcawg_expo$rel_activity[, -1])

Signature clustering

We can take a try to cluster signatures based on their relative activities.

expo_mat_rel <- pcawg_expo$rel_activity %>%
  tibble::column_to_rownames("sample") %>%
  as.matrix() %>%
  t()

hc <- hclust(dist(expo_mat_rel))

factoextra::fviz_dend(hc, k = 3, rect = T)

CNS and SBS

CNS and DBS

CNS and ID

Copy number signatures

Cancer subtyping

Prognosis